IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
Image Deblurring Using Derivative Compressed
Sensing for Optical Imaging Application
Mohammad Rostami, Student Member, IEEE, Oleg Michailovich, Member, IEEE, and Zhou Wang, Member, IEEE
o
00
<
O
>
o
o>
in
l>
o
X
Abstract — Reconstruction of multidimensional signals from the
samples of their partial derivatives is known to be a standard
problem in inverse theory. Such and similar problems routinely
arise in numerous areas of applied sciences, including optical
imaging, laser interferometry, computer vision, remote sensing
and control. Though being ill-posed in nature, the above problem
can be solved in a unique and stable manner, provided proper
regularization and relevant boundary conditions. In this paper,
however, a more challenging setup is addressed, in which one
has to recover an image of interest from its noisy and blurry
version, while the only information available about the imaging
system at hand is the amplitude of the generalized pupil function
(GPF) along with partial observations of the gradient of GPF's
phase. In this case, the phase-related information is collected
using a simplified version of the Shack-Hartmann interferometer,
followed by recovering the entire phase by means of derivative
compressed sensing. Subsequently, the estimated phase can be
combined with the amplitude of the GPF to produce an estimate
of the point spread function (PSF), whose knowledge is essential
for subsequent image deconvolution. In summary, the principal
contribution of this work is twofold. First, we demonstrate how to
simplify the construction of the Shack-Hartmann interferometer
so as to make it less expensive and hence more accessible. Second,
it is shown by means of numerical experiments that the above
simplification and its associated solution scheme produce image
reconstructions of the quality comparable to those obtained using
dense sampling of the GPF phase.
Index Terms — Deconvolution, inverse problem, derivative com-
pressive sampling, and Shack-Hartmann interferometer.
I. Introduction
THE necessity to recover digital images from their dis-
torted and noisy observations arises in a multitude of
practical applications, with some specific examples including
image denoising, super-resolution, image restoration, and wa-
termarking, just to name a few |[l)-|j4J. In such cases, it is
standard to assume that the observed image v is formed by
subjecting the original image u to convolution with a point
spread functiorT] (PSF) i, followed by contamination by white
Gaussian noise (WGN) v. Thus, formally.
i-^ u + V.
(1)
While u and v can be regarded as general members of the
signal space L2(ri) of real-valued functions on 51 C M?, the
PSF i is normally a much smoother function, with effectively
band-limited spectrum. As a result, the convolution with i has
a destructive effect on the informational content of u, in which
All the authors are with the Department of Electrical and Computer
Engineering, University of Waterloo, N2L 3G1 Ontario.
' Note that, in optical imaging, this function is also referred to as an impulse
transfer function |5 1.
case V typically has a substantially reduced set of features with
respect to u. This makes the problem of reconstruction of u
from V a problem of significant practical importance 161.
Reconstruction of the original image u from v can be carried
out within the framework of image deconvolution, which is a
specific instance of a more general class of inverse problems
JT). Most of such methods are Bayesian in nature, in which
case the information lost in the process of convolution with i
is recovered by requiring the optimal solution to reside within
a predefined functional class |8 1, |9|. Thus, for example, in the
case when u is known to be an image of bounded variation, the
above regularization leads to the famous Rudin-Osher-Fatemi
reconstruction scheme, in which u is estimated as a solution
to the following problem 1 10 1, 1 1 1 1
arg mm
2"
\l+a
|Vu| dxdy >
(2)
where a > is the regularization parameter. It should be noted
that, if the PSF obeys J i dxdy ^ 0, the problem O) is strictly
convex and therefore admits a unique minimizer, which can
be computed by a spectrum of available algorithms pO) , pTj .
A particularly non-trivial version of deconvolution is com-
monly referred to as blind. In this case, the original image
u is to be estimated without the knowledge of the PSF I?).
In this paper, however, we follow the philosophy of hybrid
deconvolution |12|, which takes advantage of any partial
information on the PSF to improve the image reconstruction.
Thus, in the algorithm described in this paper, the original
image u will be recovered from v and some partial information
on i.
Optical (and, in particular, turbulent) imaging is unarguably
the field of applied sciences from which the notion of decon-
volution has originally emanated p3j-p7j. In short-exposure
imaging, however, computational methods of image restoration
are still superseded by adaptive optics. As recently as a
decade ago, the use of adaptive optics would have been
considered as the only practical option. Nowadays, however,
with the advent of distributed cluster computing and GPU-
based image processing, it seems to be time to revisit the cost-
to-performance characteristics of the existing tools of adaptive
optics. Thus, in this work, our focus is on a specific tool of
adaptive optics, known as the Shack-Hartmann interferometer
1 16 1, p7| . Instead of completely excluding the interferometer
from our measurement system, we propose to modify its
construction through reducing the number of its local wave-
front lenses. Although the advantages of such a simplification
are immediate to see, its main shortcoming is obvious as
well: the smaller the number of lenses is, the stronger is
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
the effect of undersampling and aliasing. Accordingly, to
overcome this problem, we propose to augment the modified
Shack-Hartmann interferometer by subjecting its output to the
derivative compressed sensing (DCS) algorithm of 1 18 1. As it
will be shown later in the paper, the PSF i is determined by
a generalized pupil function P, which can be expressed in a
polar form as P = A e^"^. While the amplitude A can be mea-
sured via calibration or computed as a function of the aperture
geometry, the phase is often influenced by environmental
effects and hence it needs to be recovered from observations.
It will be shown below that DCS is particularly well suited
for reconstruction of </) from incomplete measurements of
its partial differences. Such an estimate can be subsequently
combined with A to yield an estimate of the PSF i, which
can in turn be used by a deconvolution algorithm. Thus, the
proposed method for estimation of the PSF and subsequent
deconvolution of u can be regarded as a hybrid deconvolution
technique, which comes to simplify the design and complexity
of adaptive optics on the one hand, and to make the process
of reconstruction of optical images as automatic as possible,
on the other hand.
The rest of the paper is organized as follows. Section II
summarizes basic technical preliminaries. In Section III, we
describe the SH interferometer as well as phase measurements
in optical imaging. In Section IV, we explain DCS and our new
approach to solve it. In Section V we describe deconvolution
process to recover the original image. Experimental results are
presented in Section VI, while Section VII finahzes the paper
with a discussion and main conclusions.
II. Technical Preliminaries
In short exposure imaging, due to aberrations in the imaging
system induced by, e.g., atmospheric turbulence, the impulse
response of an optical imaging system is often unknown p9) .
In order to better understand the setup under consideration, we
first note that, in optical imaging, the PSF i is obtained from
an amplitude spread function (ASF) h as i :— |ft,p. The ASF,
in turn, is defined in terms of the generalized pupil function
(GPF) P(x, y) as given by 1,20}
K^.ii)^
1
A.
P{x,y)e
rix^+yri)
dxdy, (3)
where Zi is the focal distance and A^, is the optical wavelength.
Being a complex-valued quantity, P{x, y) can be represented
in terms of its amplitude A{x, y) and phase (/)(x, y) as
P{x,y):^A{x,y)e^'f''-'-''y\
(4)
Here, the GPF amplitude A{x,y) (which is sometimes simply
referred to as the aperture function) is normally a function
of the aperture geometry. Thus, for instance, in the case of a
circular aperture, A{x, y) can be defined as ||19|
A{r) =
ifr<f
otherwise
(5)
where D denotes the pupil diameter. Thus, given (j){x,y), one
could determine h and therefore i. Unfortunately, the phase
(j){x, y) does not have an analytic expression, and it has to be
measured in practice using such tools as the Shack-Hartmann
interferometer (SHI) | fT6l .
As will be discussed later in the paper, the SHI is capable
of sensing the partial derivatives of (t>{x,y). Needless to say,
in order to minimize the effect of aliasing on the estimation
result, an accurate reconstruction of 0(a:, y) requires taking a
fairly large number of the samples of V0(x, y) |21 1. In some
applications, the number of sampling points (as defined by the
number of local wavefront lenses) reaches as many as a few
thousands. It goes without saying that reducing the number of
lenses would have a positive impact on the SHI in terms of its
cost and approachability. Alas, such a reduction is impossible
without undersampling, which tends to have formidable effect
on the overall quality of phase estimation.
In this paper, to minimize the effect of undersampling, we
exploit DCS |18|. As opposed to the classical compressed
sensing (CCS) |22| , in addition to the sparsifing constraints,
DCS also uses constraints which are intrinsic in the definition
of partial derivatives. Using these additional constraints -
which are called the cross-derivative constraints - allows sub-
stantially improving the quality of reconstruction of (t)[x,y),
as compared to the case of CCS-based estimation.
III. Shack-Hartmann Interferometer (SHI)
As it was mentioned earlier, the SHI is typically used to
measure the gradient V0(x, y) of the GPF phase 0(a;, y), from
which the values of the latter can be subsequently estimated.
To this end, the unknown phase (f){x, y) is assumed to be
expandable in terms of some basis functions {Zk}'^^^, viz.
^ ' oo
(l^{x,y)^^akZk{x,y), (6)
fc=0
where the representation coefficients {ak}kLo are assumed
to be unique and stably computable. Note that, in this case,
the datum of {flfej^Q uniquely identifies (f){x,y), while the
coefficients {a/d^g can be estimated due to the linearity of
(|6]l which suggests
V(?!)(x, y) = ^ak ^Zk{x, y),
(7)
fe=0
The most frequent choice of {^fe}^o ^'^ adaptive optics is
the Zernike polynomials (aka Zemike functions) [20 1. These
polynomials constitute an orthonormal basis in the space of
square-integrable functions defined over the unit disk in M^.
Zernike polynomials can be subdivided in two subsets of the
even Z™ and odd Z,7™ Zernike polynomials which have very
convenient analytical definitions as given by
Z::^ip,^)^R^{p)cos{mip)
(8)
(9)
where m and n are nonnegative integers with n > m, <
If < 271 is the azimuthal angle, and < p < 1 is the radial
distance. The radial polynomials _R™ are defined as
(n-m)/2
Kip)
(-!)'=(« -fc)!
^ A:!((n + m)/2-/s)!((n-m)/2-fc)!
11-2 k
(10)
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
Fig. 1. An example of a 10 X 10 SHI array on a circular aperture. The
shading indicates those blocks (i.e., lenses) which are rendered idle.
Note that, since the Zernike polynomials above are defined
using polar coordinates, it makes sense to re-express the
phase (p and its gradient in the polar coordinate system as
well. (Technically, this would amount to replacing x and y
in (|6]l-(j7]i by p and ip, respectively.) Moreover, due to the
property of the Zernike polynomials to be an orthonormal
basis, the representation coefficients {a^l^Q in in ((6])-(|7]i can
be computed by orthogonal projection, namely
a-k
(t){p,(p) Zk{p,(p) pdpdip
(11)
In practice, however, (j){p, (p) is unknown and therefore the
coefficients {flfej^o n^ed to be estimated by other means.
Thus, in the case of the SHI, the coefficients can be estimated
from a finite set of discrete measurements of V0(p, (p).
The main function of the SHI is to acquire discrete measure-
ments of V0 by means of linearization. The linearization takes
advantage of subdividing a (circular) aperture into rectangular
blocks with their sides formed by a uniform rectangular lattice.
An example of such a subdivision is shown in Fig. [Tlfor the
case of a 10 X 10 lattice grid. Subsequently, it is assumed
that the grid is sufficiently fine to approximate a restriction of
the phase <j) to each of the above blocks by a linear function.
This results in a piecewise linear approximation of 0, whose
accuracy improves asymptotically when the lattice size goes
to infinit>|^ Formally, let n := {{x,y) e M^ | x^ + y^ < jj2y
be a circular aperture of radius D and 5 = {(x, y) e M^ |
max{|a;|, \y\} < D} be a square subset of R-^ such that fi C 5.
Then, for each polar coordinate (p, ip) £ ft and an N x N grid
of square blocks of size 2D/N x 2D/N, the phase (j) can be
expressed as
4>(x,y) ~ ax + by + c^ (12)
for all (x, y) in a neighbourhood of (p cos (p,p sin ip). The
approximation in ([T2| suggests that
V0(x,y)« [a, 6]^
(13)
^More rigorously, one can show that, as long as </> is uniformly continuous,
its piecewise linear approximation converges uniformly, as the grid size goes
to infinity.
LoisAnij C CD Clip
Shack-Haattraan Pattern
Fig. 2. Basic structure of the SHI and a resulting pattern of the focal points.
where (•)-^ denotes matrix transposition. While c in ( 12 1 can be
derived from boundary conditions, coefficients a and b should
be determined via direct measurements. To this end, the SHI
is endowed with an array of small focusing lenses, which are
supported over each of the square blocks of the discrete grid. In
the absence of phase aberrations, the focal points of the lenses
are spatially identified and registered using a high-resolution
CCD detector, whose imaging plane is aligned with the plane
of the focal points. Then, when the wavefront gets distorted
as a result of, e.g., atmospheric turbulence, the focal points
are "pushed" towards new spatial positions, which can also be
pinpointed by the same detector. The resulting displacements
can therefore be measured and subsequently related to the
values of V0 at corresponding points.
To explain how the above procedure can be performed,
additional notations are in order. Let ft^ denote a finite set
of spatial coordinates defined according to
n.
Xd
yd
\{xd,yd)
e n
-D+^(i + ^], t^O,l,...,N~l (14)
and xl + yj<D^y
Note that the set fid can be thought of as a set of the
spatial coordinates of the geometric centres of the wavefront
lenses, restricted to the domain of aperture O. Under some
reasonable assumptions, one can then show 1 ,23) that the focal
displacement A{x,y) — [Ax{x, y), Ay{x,y)Y' measured at
some {x, y) G fid is related to the value of V(/)(x, y) as given
by
V0(a;, y) « -A^{x, y), V(a;, y) G ^d, (15)
where / is the focal distance of the wavefront lenses. Such a
measurement setup is depicted in Fig.|2]along with an example
of measured focal points.
Now, given a total of M :— #51^ measurements of V0
over rijj, one can try to recover a useful approximation of (j)
over the whole Vl in the form of a projection of 4> on the
linear subspace spanned by all Zernike polynomials up to the
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
order L inclusive. In this case, it is possible to estimate the
representation coefficients {ak\k=a ^^^ solution of
conditions of compressed sensing far superseded those of the
subject to appropriate boundary conditions. It is worthwhile
noting that (16 1 can be rewritten in a vector-matrix form as
IZa-dl
2
2j
s.t. a h 0,
(17)
where Z is a 2Af x i + 1 matrix composed of the values
of the partial derivatives of the Zernike polynomials, d is
a measurement (column) vector of length 2M, and a —
[oq, fli, . . . , ai]^ is a vector of the representation coefficients
of (j). The constraint a )^ in (17 1 is supposed to further
regularize the solution by forcing a to be a member of some
convex cone as well. Thus, for example, if the mean value of
(j) can be assumed to be equal to zero, the solution to ( 17 1 can
be computed as
a=Z#d, (18)
where Z* denotes the pseudo-inverse of Z, whose definition
is unique and stable as long as the row-rank of Z is greater or
equal to L + 1 (hence suggesting that 2M > L + 1). Having
estimated a, the phase can be approximated as
<P{P,'~P) ~ ^akZk(p,(p).
(19)
k=0
A higher accuracy of phase estimation requires using higher-
order Zernike polynomials, which in turn necessitates a pro-
portional increase in the number of wavefront lenses. More-
over, as required by the linearization procedure in the SHI,
the lenses have to be of a relatively small sizes (sometimes,
on the order of a few microns), which may lead to the use
of a few thousand lenses per one interferometer. Accordingly,
to simplify the construction and to reduce the cost of SHIs,
we propose to substantially reduce the number of wavefront
lenses, while compensating for the induced information loss
through the use of DCS, which is detailed next.
IV. Derivative Compressive Sampling
A. Classical Compressed Sensing
Central to signal processing is the Shannon-Nyquist theorem
p4) , which specifies conditions on which a band-limited
signal can be stably and uniquely recovered from its discrete
measurements. However, in around 2005, a different sampling
theorem was formed which, in some cases, abrogates the
fundamentals of its predecessor. This new theory, nowadays
known as compressed sensing (aka compressive sampling),
asserts that signals, which admit a sparse representation in
a predefined basis/frame, can be recovered from their discrete
measurements, whose number is proportional to the ^o-norm
of the coefficients of the sparse representation. In such a
case, the sparser the representation of the signal is, the
smaller can be the number of measurements required for signal
reconstruction. As a result, cases are numerous in which the
Shannon-Nyquist sampling |22|, |25|.
Despite its widespread success in countless applications, the
theory of compressed sensing is not entirely free of limitations.
One of such limitations stems from the necessity to use a non-
linear decoder. Indeed, while in the case of classical sampling,
the reconstruction of a time-domain signal is implemented
through linear interpolation (i.e., linear filtering), in the case
of compressed sensing, the reconstruction involves solution of
a convex optimization problem. Specifically, let us consider a
typical setup of classical compressed sensing (CCS), in which
y £ M'" represents an observed version of x e M", related
according to
y = *x, (20)
where ^ G M™^" is an observation (sampling) matrix with
n > m.
The recovery of x from y based on ( [20| is impossible to
implement in a unique and stable way, unless it is known that
X is sparse and hence has a relatively low value of ||x||o. In
such a case, if the sampling matrix ^ satisfies the restricted
isometry property (RIP) f22|, p5) with respect to a certain
class of sparse signals to which x is assumed to belong, then
CCS recovers x as a solution to | |26l , | |27l
X = argmin{||x'||i | ^x' = y} , (21)
x'
which is a convex minimization problem, which is straightfor-
ward to reformulate in terms of linear programming. Moreover,
in the case when the measurements y are error-prone, a more
robust version of CCS is to recover x as given by
X — arg mm
{llx'l
|*x'
y\\l
<^}.
(22)
where e > is a parameter controlling the size of the noise.
Moreover, it was shown in 1221, ||25l, that the estimation error
in the signal reconstructed according to ( 22 1 can be bounded
by a linear function of e. This implies robustness of the CCS
reconstruction towards the presence of measurement noise.
It should be finally noted that the optimization problem
(22 1 can be reformulated in its alternative Lagrangian form.
in which case one can find
1
X — arg mm
*x'
y||^
Allx'l
(23)
where A > is an optimal Lagrange multiplier p6) . In what
follows, it is assumed that an optimal value of A is known.
(For more details on this subject, the reader is referred to p6|
as well as to the later sections of this paper).
B. Derivative Compressed Sensing (DCS)
Let the partial derivatives of evaluated (by means of the
SHI) at the points of set ilij be column-stacked into vectors
fj; and iy of length M = i^fld- In what follows, the partial
derivatives i^ and fj, are assumed to be sparsely representable
by an orthonormal basis in M*^. Representing such a basis by
an M X M unitary matrix W, the above assumption suggests
the existence of two sparse vectors c^, and Cy such that
Wcx and Wcy provide accurate approximations of the partial
derivatives of the original phase (j) evaluated over fid-
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
Now, the simplification of the SHI proposed in the current
paper amounts to reducing the number of wavefront lenses
to a minimum. Formally, such a reduction can be described
by two n X M sub-sampling matrices '^^ and '^y, where
n < M. Specifically, let h^ '■— 'i'xix and bj, :— 'i'yiy be
incomplete (partial) observations of i^ and fj^, respectively.
Then, the noise-free counterparts of the partial derivatives can
be approximated by Wc^ and Wc* respectively, where
< = argmin |i||*,W^< - b,||2 + A,||<||i| (24)
where p^*) is a vector of Bregman variables (or, equivalently,
augmented Lagrange multipliers) and (5 > is a user-defined
parameteipl
and
argmin<; -W'i'yWc'y - hy\\l
^v\W%
Vll^yl
(25)
for some Aa;,Ay > 0. Moreover, in the case when A^^ — Xy,
the above estimates can be combined together To this end, let
c = [Cx,Cyf, b = [K,hyf, and A = diag{*^M^, *yTy} e
M2"x2M^ Then,
argmm ■^ -||Ac'
b||2 + A||c'||i
(26)
The c-update step in ( 32 1 has the format of a standard
basis pursuit de-noising (BPDN) problem f29l, which can
be solved by a variety of optimization methods [301. In the
present paper, we used the FISTA algorithm of |31) due to the
simplicity of its implementation as well as for its remarkable
convergence properties. It should be noted that the algorithm
does not require explicitly defining the matrices A and B.
Only the operations of multiplication by these matrices and
their transposes need to be known, which can be implemented
in an implicit and computationally efficient manner
Once an optimal c* is recovered, it can be used to estimate
the noise-free versions of i^ and fy as WT^c* and WTyC*,
respectively. These estimates can be subsequently passed on
to the fitting procedure of Section III to recover the values of
(j), which, in combination with a known aperture function A,
provide an estimate of the PSF i as an inverse discrete Fourier
transform of the autocorrelation of P = AeJ'^. Algorithm 1
below summarizes our method of estimation of the PSF.
where A = A, = Xy. In this form, the problem (|26f is identical Algorithm 1: PSF estimation via DCS
to (23 I and hence it can be solved by a variety of optimization
algorithms p6), p7).
The DCS algorithm augments CCS by subjecting the mini-
mization in (|26| to an additional constraint which stems from
the fact that |18|
dx dy dy dx '
which is valid for any two times continuously differentiable
4i{x^ y). Thus, in particular, the constraint implies the existence
of two partial differences matrices D^ and Dy which obey
DJy = Dyix. (28)
Consequently, if T^. and Ty are the matrices satisfying T^c —
1) Data: h^, by, and A >
2) Initialization: For a given transform matrix W and
matrices/operators '^x, ^y, D^, Dy, T^ and Ty, preset
the procedures of multiplication by A, A^ , B and B^ .
3) Phase recovery: Starting with an arbitrary c^^^ and
p'-''-' = 0, iterate ( 32 1 until convergence to result in an
Cx and TyC = Cy, respectively, then (31 1 can be re-expressed
in terms of c as
(29)
or
^y-^x^ — J-^x-^y^
Bc = 0,
(30)
where B :— DyTx — DxTy. Thus, DCS solves the constrained
minimization problem as given by
c* = argmin |i||Ac' ~h\\l + A||c'||i| , (31)
s.t. Be' =
optimal c*. Use the estimated (full) partial derivatives
WTxC* and WTyC* to recover the values of (p over H..
4) PSF estimation: Using a known aperture function A,
compute the inverse Fourier transform of P = A eJ^ to
result in a corresponding ASF h. Estimate the PSF i as
^^\h\\
The estimated PSF can be used to recover the original image
u from V through the process of deconvolution as explained
in the section that follows.
V. Deconvolution
The acquisition model ([TJ can be rewritten in an equivalent
operator form as given by
V ~ H{u} + V,
(33)
A solution to ( 3 1 1 can be found, for instance, by means of
the Bregman algorithm L28J, in which case c* is approximated
by a stationary point of the sequence of iterations produced
by
^(t+i) ^
arg miiic'
(211^'- - 1^112-
\Ac'-h\\
+A||c'||i
I Be'
.(*)
IB}
(32)
where Ti, denote the operator of convolution with the estimated
PSF i. Note that, in this case, the noise term v accounts for
both measurement noise as well as the inaccuracies related to
estimation error in i.
The deconvolution problem of finding a useful approxima-
tion of u given its distorted measurement v can be addressed in
many way, using a multitude of different techniques pTj-p3).
In this work, we use the ROF model and recover a regularized
^In this work, we use 5 = 0.5.
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
approximation as a solution of
u* = argniin <^ -\\n{u} - v\\l + 7 ||u|JTy
u I Z
(34)
where ||u||tv = J J\Wu\dxdy denotes the total variation
(TV) semi-norm of u.
One computationally efficient way to solve ( (34] i is to
substitute a direct minimization of the cost function in (|34])
by recursively minimizing a sequence of its local quadratic
majorizers |31|. In this case, the optimal solution u* can
be approximated by the stationary point of a sequence of
intermediate solutions produced by
w
(«)
,(*)
M'H*{w--H{w(*'}}
,(t+i) ^ argmin„ {^\\u ~ u.(*)||2 + 7 \\u\\tv} ,
(35)
where H* is the adjoint of H and fi is chosen to satisfy /i >
IIH*?^!!. In this paper, the TV denoising at the second step
of ([35| has been performed using the fixed-point algorithm
of Chambolle fTTl. The convergence of ( |35] l can be further
improved by using the same FISTA algorithm of pT| . The
resulting procedure is summarized below in Algorithm 2.
Algorithm 2: TV deconvolution using FISTA
1) Initialize: Select an initial value m'^^; set y'"^ = u^^^^
and t(o) = 1
2) Repeat until convergence:
w
(t+i)
y^*'> + ^lH* {v - H{y^*^}}
argmin„{i||M-u;(*)|||+7||u|JTy}
(*+i) = 0.5 (1 + v'1 + 4(t(*))2)
,(*+!)
In summary. Algorithms 1 and 2 represent the essence of
the proposed algorithm for hybrid deconvolution of optical
images. The next section provides experimental results which
further support the value and applicability of the proposed
methodology.
VI. Results
To demonstrate the viability of the proposed approach, its
performance has been compared against reference methods.
The first reference method used a dense sampling of the phase
(as it would have been the case with a regular SHI), thereby
eliminating the need for a CS-based phase reconstruction. The
resulting method is referred below to as the dense sampling
(DS) approach. Second, to assess the importance of incorpora-
tion of the cross-derivative constraints, we have used both CCS
and DCS for phase recovery. In what follows, comparative
results for phase estimation and subsequent deconvolution are
provided for all the above methods.
A. Phase recovery
To assess the performance of the proposed and reference
methods under controllable conditions, simulation data was
used. The random nature of atmospheric turbulence necessi-
tated the use of statistical methods to model its effect on a
(a)
(b)
(c)
Fig. 3. An example of a simulated phase (j> (subplot (a)) along with its partial
derivatives w.r.t. x (subplot (b)) and y (subplot (c)).
wavefront propagation. Specifically, in this paper, the effect of
atmospheric turbulence has been described using the modified
Von Karman PSF model |34|. A typical example of a GPF
phase (p is shown in subplot (a) of Fig. [5] In the shown case,
the size of the phase screen was set to be equal to 10 x 10 cm,
while the sampling was performed over a 128 x 128 uniform
lattice (which would have corresponded to the use of 16384
lenses of a SHI). The partial derivatives d(j)/dx and d(j)/dy
are shown in subplots (b) and (c) of Fig. [3] respectively.
In the present paper, the subsampling matrices '^^ and ^y
were obtained from an identity matrix / through a random
subsampling of its rows to result in a required compression
ratio r (to be specified below). To sparsely represent the partial
derivatives of (j)^ ^ was defined to correspond to a four-
level orthogonal wavelet transform using the nearly symmetric
wavelets of Daubechies with five vanishing moments |35J.
To demonstrate the value of using the cross derivative con-
straint for phase reconstruction, the CCS and DCS algorithms
have been compared in terms of the mean square errors (MSE)
of their corresponding phase estimates. The results of this
comparison are summarized in Fig. |4]for different compression
ratios (or, equivalently, (sub)sampling densities) and SNR = 40
dB.
As expected, one can see that DCS results in lower values
of MSE as compared to CCS, which implies a higher accuracy
of phase reconstruction. Moreover, the difference in the per-
formances of CCS and DCS appears to be more significant for
lower sampling rates, while both algorithms tend to perform
similarly when the sampling density approaches the DS case.
Specifically, when the sampling density is equal to r = 0.3,
DCS results in a ten times smaller MSE as compared to
the case of CCS, whereas both algorithms have comparable
performance for r — 0.83. This result characterizes DCS as
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
Q H n
— []
— ei
Fig. 4. MSB of phase reconstruction obtained with different methods as a
function of r. Here, the dashed and solid lines correspond to CCS and DCS,
respectively, and SNR is equal to 40 dB.
a better performer than CCS in the case of relatively low
(sub)sampling rates.
Some typical phase reconstruction results are shown in
Fig.|5] whose left and right subplots depict the phase estimates
obtained using the CCS and DCS algorithms, respectively, for
the case of r = 0.5. To visualize the differences more clearly,
error maps for both methods are shown in subplot (c) and
(d)(note for better visualization error map values are multiplied
by factor of 10.). A close comparison with the original phase
(as shown in subplot (a) of Fig. [3]) reveals that DCS provides
a more accurate recovery of the original (f), which further
supports the value of using the cross-derivative constraints. In
fact, exploiting these constraints effectively amounts to using
additional measurements, which are ignored in the case of
CCS.
To investigate the robustness of the compared algorithms
towards the influence of additive noises, their performances
have been compared for a range of SNR values. The results
of this comparison are summarized in Fig. [6] Since the cross-
derivative constraints exploited by DCS effectively restrict
the feasibility region for an optimal solution, the algorithm
exhibits a substantially better robustness to the additive noise
as compared to the case of CCS. This fact represents another
beneficial outcome of incorporating the cross-derivative con-
straints in the process of phase recovery.
It should be taken into account that, although the shape of
does not change the energy of the PSF i, it plays a crucial role
in determining its spatial behaviour In the section that follows,
it will be shown that even small inaccuracies in reconstruction
of (j) could be translated into dramatic difference in the quality
of image deconvolution.
B. Image deconvolution
As a next step, the phase estimates obtained using the CCS-
and DCS-based methods for r — 0.5 were combined with the
aperture function A to result in their respective estimates of the
PSF i. These estimates were subsequently used to deconvolve
(c)
(d)
Fig. 5. (Subplot (a)) Phase reconstructed obtained by means of CCS for
SNR = 40 dB and r = 0.5; (Subplot (b)) Phase reconstructed obtained by
means of DCS for the same values of SNR and r.; (Subplot (c) and (d))
Corresponding error maps for CCS and DCS.
Fig. 6. MSE of phase reconstruction obtained with different methods as a
function of SNR. Here, the dashed and solid lines con'espond to CCS and
DCS, respectively, and r = 0.5.
a number of test images such as "Satellite", "Saturn", "Moon"
and "Galaxy". All the test images were blurred with an original
PSF, followed by their contamination with additive Gaussian
noise of different levels. As an example. Fig. p\ shows the
"Satellite" image (subplot (a)) along with its blurred and noisy
version (subplot (b)).
Using the PSF estimates, the deconvolution was carried out
using the method detailed in 1 1 1 1. For the sake of comparison,
the deconvolution was also performed using the PSF recovered
from dense sampling (DS) of (p. Note that this reconstruction
is expected to have the best accuracy, since it neither involves
undersampling nor requires a CS-based phase estimation.
All the deconvolved images have been compared with their
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
(a)
(b)
Fig. 7. Satellite image (subplot (a)) and its blurred and noisy version (subplot
(b)).
(a)
(b)
Fig. 8. (Subplot (a)) Image estimate obtained with the CCS-based method
for phase recovery (SSIM = 0.917); (Subplot (b)) Image estimate obtained
with the DCS-based method for phase recovery (SSIM = 0.781).
has been shown to yield phase estimates of substantially better
quality as compared to the case of CCS.
In this paper, our main focus has been on simplifying
the structure of the SHI through reducing the number of
its wavefront lenses, while compensating for the effect of
undersampling by using the theory of CS augmented by the
cross-derivative constraint. The solution was computed using
the Bregman algorithm, which provides a computationally ef-
ficient framework to carry out the constrained phase recovery.
Moreover, the resulting phase estimates were used to recover
their associated PSF, which was subsequently used for image
deconvolution. It was shown that the DCS-based estimation of
(j) with r = 0.3 results in image reconstructions of the quality
comparable to that of DS, while substantially outperforming
the results obtained with CCS.
While the proposed method offers a practical solution to
the problem of phase estimation in adaptive optics, some
interesting questions about the theoretical aspects of DCS still
lay open. In particular, the question of theoretical performance
of CS in the presence of side information on the source signal
needs to be addressed through future research.
Acknowledgment
This work was supported in part by the Natural Sciences
and Engineering Research Council of Canada and in part by
Ontario Early Researcher Award program, which are gratefully
acknowledged. The authors would also like to acknowledge
Sudipto Dolui for his helpful comments as well as for provid-
ing deconvolution codes.
original counterparts in terms of PSNR as well as of the
structural similarity index (SSIM) of |36|, , which is believed
to be a better indicator of perceptual image quality |37|. The
resulting values of the comparison metrics are summarized in
Table 1, while Fig.lSlshows the deconvolution results produced
by the CCS- and DCS-based methods.
The above results demonstrate the importance of accurate
phase recovery, where even a relatively small phase error can
have a dramatic effect on the quality of image deconvolution.
Under such conditions, the proposed method produces image
reconstructions of a superior quality as compared to the case of
CCS. Moreover, comparing the results of Table 1, one can see
that DS only slightly outperforms DCS in terms of PSNR and
SSIM, while in many practical cases, the difference between
the performances of these methods are hard to detect visually.
VII. Discussion and conclusions
In the present paper, the applicability of DCS to the problem
of reconstruction of optical images has been demonstrated. It
was shown that, in the presence of atmospheric turbulence,
the phase <j) of the GPF P = AeJ'^ is a random function,
which needs to be measured using adaptive optics. To simplify
the complexity of the latter, a CS-based approach has been
proposed. As opposed to CCS, however, the proposed method
performs phase reconstruction subject to an additional con-
straint, which stems from the property of V0 to be a potential
field. The resulting algorithm (referred to as the DCS method)
References
[1] J. Yang, J. Wright, T. S. Huang, and Y. Ma, "linage super-resolution
via sparse representation," IEEE Transactions on Image Processing, vol.
19, pp. 2861-2873, 2010.
[2] M. Elad and M. Aharon, "Image denoising via sparse and redundant
representations over learned dictionaries," IEEE Transactions on Image
Processing, vol. 17. pp. 3736-3745, 2006.
[3] IJ. Mairal, G. Sapiro, and M. Elad, "Learning multiscale sparse
representations for image and video restoration," Multiscale Modeling
and Simulation, vol. 7, pp. 214-241, 2008.
[4] R. T. Paul, "Review of robust video watermarking techniques," IJCA
Special Issue on Computational Science, , no. 3, pp. 90-95, 2011.
[5] G. D. Boreman, Modulation Transfer Function in Optical and Electro-
Optical Systems, SPIE Optical Engineering Press, Bellingham, Wash-
ington, 2001.
[6] R. T. Paul, "Blind deconvolution via cumulant extrema," IEEE Signal
Processing Magazine, , no. 3, pp. 2A-A2, 1996.
[7] D. Kundur and D. Hatzinakos, "Blind image deconvolution," IEEE
Signal Processing Magazine, , no. 3, pp. 43-64, 1996.
[8] S. Geman and D. Geman, "Stochastic relaxation, gibbs distribution and
the bayesian restoration of images," IEEE Trans. Pattern Analysis and
Machine Intelligence, vol. PAMI-6, pp. 721-741, 1984.
[9] V. Ton'e T. Poggio and C. Koch, "Computational vision and regulariza-
tion theory," Nature, vol. 317, pp. 314-319, 1985.
[10] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based
noise removal algorithms," Phys. D, vol. 60, pp. 259-268, November
1992.
[11] A. Chambolle, "An algorithm for total variation minimization and
applications," Journal of Mathematical Imaging and Vision, vol. 20,
pp. 89-97, 2004.
[12] O. Michailovich and A. Tannenbaum, "Blind deconvolution of medical
ultrasound images: Parametric inverse filtering approach," IEEE Trans-
actions on Image Processing, vol. 16, no. 12, pp. 3005-3019, December
2007.
[13] W. H. Richardson, "Bayesian-based iterative method of image restora-
tion," / Opt. Soc. Am. A, vol. 62, no. 1, pp. 55-59, 1972.
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. XX, NO. XX, JANUARY 2012
TABLE I
SSIM AND PSNR COMPARISONS OF PHASE RECOVERY RESULTS
Image
Satellite
Saturn
Moon
Galaxy
Noise std
10-^
0.001
0.003
0.005
10--''
0.001
0.003
0.005 1 10-^
0.001
0.003
0.005
10-'''
0.001
0.003
0.005
PSNR comparison (in
dB)
Blurred
14.06
14.06
14.06
14.05
17.78
17.78
17.78
17.78
19.98
19.97
19.97
19.97
18.79
18.79
18.78
18.78
DS
27.97
27.75
25.97
22.43
31.49
31.08
28.50
23.89
25.06
25.04
24.83
21.76
23.58
23.60
23.38
20.93
CS
17.06
16.93
16.54
15.63
23.42
23.38
22.80
20.55
22.36
22.38
22.30
19.73
21.16
21.12
20.64
18.46
DCS
27.42
27.22
25.56
22.22
31.02
30.65
28.30
23.72
25.00
24.99
24.78
21.73
23.52
23.54
23.32
20.86
SSIM comparison
Blurred
0.200
0.200
0.199
0.197
0.226
0.226
0.226
0.175
0.512
0.512
0.509
0.504
0.257
0.257
0.257
0.254
DS
0.730
0.720
0.554
0.269
0.688
0.660
0.506
0.228
0.645
0.642
0.607
0.552
0.493
0.495
0.501
0.397
CS
0.349
0.344
0.306
0.206
0.424
0.416
0.348
0.212
0.539
0.538
0.493
0.488
0.348
0.347
0.326
0.224
DCS
0.674
0.667
0.519
0.263
0.656
0.641
0.483
0.223
0.643
0.640
0.604
0.549
0.490
0.491
0.501
0.393
[14] L. B. Lucy, "An iterative technique for the rectification of observed
distributions," Astron. /, vol. 79, no. 6, pp. 745-754, 1974.
[15] P. A. Jansson, Deconvolution of Images and Spectra, Opt. Eng. 36,
3224, 1997.
[16] D. Dayton, B. Pierson, B. Spielbusch, and J. Gonglewski, "Atmospheric
structure function measurements with a shack-hartmann wave-front
sensor," Journal of Mathematical Imaging and Vision, vol. 20, pp. 89-
97, 2004.
[17] R. G. Lane and M. Tallon, "Wave-front reconstruction using a shack-
hartmann sensor," Applied Optics, vol. 31, pp. 6902-6908, 1992.
[18] M. Hosseini and O. Michailovich, "Derivative compressive sampling
with application to phase unwrapping," in Proceedings of EUSIPCO,
Glasgow, UK, August 2009.
[19] M. C. Roggemann and B. M. Welsh, Imaging Through Turbulence, CRC
Press, Boca Raton, FL, 1996.
[20] D. L. Fried, "Statistics of a geometric representation of wavefront
distortion," J. Opt. Soc. Am., vol. 55, pp. 1427-1431, 1965.
[21] O. Michailovich and A. Tannenbaum, "A fast approximation of smooth
functions from samples of partial derivatives with application to phase
unwrapping," Signal Processing, vol. 88, pp. 358-374, 2008.
[22] Y. Tsaig and D. L. Donoho, "Compressed sensing," IEEE Trans. Inform.
Theory, vol. 52, pp. 1289-1306, 2006.
[23] J. Primot, G. Rousset, and J. C. Fontanella, "Deconvolution from wave-
front sensing: a new technique for compensating turbulence-degraded
images," / Opt. Soc. Am., vol. 7, pp. 1598-1608, 1990.
[24] C. E. Shannon, "Communication in the presence of noise," Proceedings
of the IRE, vol. 37, no. 1, pp. 10-21, 1949.
[25] E. J. Candes, J. Romberg, and T. Tao, "Robust uncertainty principles:
exact signal reconstruction from highly incomplete frequency informa-
tion," IEEE Trans, hifo. Theory; vol. 52, no. 2, pp. 489-509, 2006.
[26] D. L. Donoho and Y. Tsaig, "Fast solution of ii-norm minimization
problems when the solution may be sparse," 2006.
[27] W. Yin, S. Osher, D. Goldfarb, and J. Darbon, "Bregman iterative al-
gorithms for ^ J -minimization with applications to compressed sensing,"
SIAM J. Imaging Sciences, vol. 1, no. I, pp. 143-168, 2008.
[28] J. Cai, S. Osher, and Z. Shen, "Split bregman methods and frame based
image restoration," Multiscale Modeling & Simulation, vol. 8, no. 2, pp.
337-369, 2009.
[29] S. S. Chen and D. L. Donoho, "Atomic decomposition by basis pursuit,"
Siam Journal on Scientific Computing, vol. 20, 1998.
[30] I. Daubechies, M. Defrise, and C. D. Mol, "An iterative thresholding
algorithm for linear inverse problems with a sparsity constraint," Com-
munications on Pure and Applied Mathematics, vol. 57, pp. 1413-1457,
2004.
[31] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding
algorithm for linear inverse problems," SIAM Journal on Imaging
Sciences, vol. 2, pp. 183-202, 2009.
[32] A. Savitzky and M. J. E. Golay, "Smoothing and differentiation of
data by simplified least squares procedures," Anal. Chem., vol. 36, pp.
1627-1639, 1964.
[33] A. N. Tikhonov and V. Y. Arsenin, "Solutions of ill-posed problems,"
1977.
[34] J. D. Schmidt, Numerical .simulation of optical wave propagation with
examples in MATLAB, SPIE, Washington, 2010.
[35] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Reg. Conf. Series
in Applied Math. SIAM, 1992.
[36] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. R Simoncelli, "Image
quality assessment: From error visibility to structural similarity," IEEE
Trans. Image Processing, vol. 13, no. 4, pp. 600-612, 2004.
[37] Z. Wang and A. C. Bovik, "Mean squared error: love it or leave it? - a
new look at signal fidelity measures," IEEE Signal Processing Magazine,
vol. 26, no. 1, pp. 98-117, 2009.